﻿Imports Assimp
Imports NahRT
Imports OpenTK
Imports OpenTK.Graphics.OpenGL4


Public Class SceneManager
    Public Const MaxTriangles As Integer = 16384
    Public Const Stride As Integer = 16384

    Private scene As Scene
    Public dataTexture As Integer REM 三角形数据
    Private triangleCount As Integer REM 三角形数量
    Private sceneData As IWriteSceneData
    Public ReadOnly Property TriangleData As SeparateSceneData
        Get
            Return sceneData
        End Get
    End Property

    Public Sub New(s As Scene)
        scene = s
        sceneData = New SeparateSceneData(GetTriangleCount())
    End Sub

    Public texDic As New Dictionary(Of Integer, Integer)
    Public Sub MakeSceneTexture()
        Dim texcount As Integer = 0
        texDic.Clear()
        For i = 0 To scene.MaterialCount - 1
            Dim slots = scene.Materials(i).GetMaterialTextures(TextureType.Diffuse)
            If slots IsNot Nothing Then
                Dim texID = GL.GenTexture()
                texDic.Add(i, texID)
                GL.ActiveTexture(TextureUnit.Texture0 + 5 + i) REM 保留5个内置纹理
                Dim path As String = scene.Materials(i).TextureDiffuse.FilePath
                If path IsNot Nothing Then
                    LoadTexture(path, texID)
                End If
                texcount += 1
            End If
        Next
        If texcount > 10 Then
            Throw New Exception("场景最大支持10张纹理")
        End If

        SetupSamplers()
    End Sub

    Public texturehandles() As ULong
    Public Sub MakeSceneTextureHandle()
        texDic.Clear()
        texturehandles = New ULong(scene.MaterialCount - 1) {}
        For i = 0 To scene.MaterialCount - 1
            Dim slots = scene.Materials(i).GetMaterialTextures(TextureType.Diffuse)
            Dim texID = 0
            texID = GL.GenTexture()
            texDic.Add(i, texID)
            If slots IsNot Nothing Then
                Dim path As String = scene.Materials(i).TextureDiffuse.FilePath
                If path IsNot Nothing Then
                    LoadTexture(path, texID)
                End If
            End If
            Dim hTex As ULong = GL.Arb.GetTextureHandle(texID)
            GL.Arb.MakeTextureHandleResident(hTex)
            Dim b = GL.Arb.IsTextureHandleResident(hTex)
            texturehandles(i) = hTex
        Next
    End Sub

    Private Sub SetupSamplers()
        Dim program As Integer = ShaderManager.Active("RayTracing")
        For i = 0 To 30 - 1
            Dim SamplerX As Integer = GL.GetUniformLocation(program, "Samplers[" + i.ToString() + "]")
            GL.Uniform1(SamplerX, i + 1) REM 从TextureUnit.Texture1 开始
        Next
    End Sub

    Public cameraPos As Vector4 = New Vector4(0, 0, 0, 1)
    Public mtxModel As Matrix4 = Matrix4.Identity
    Public mtxView As Matrix4 = Matrix4.Identity
    Public mtxProjection As Matrix4

    Public Sub UpdateUniform(program As Integer)
        GL.UseProgram(program) REM 注意在设置Uniform 前要先UseProgram

        Dim uMVP As Integer = GL.GetUniformLocation(program, "MVP")
        GL.UniformMatrix4(uMVP, False, mtxModel * mtxView * mtxProjection)

        'Dim glposition = Vector4.Transform(New Vector4(0, 0, 4, 1), mtxModel * mtxView * mtxProjection)

        Dim uTriangleCount As Integer = GL.GetUniformLocation(program, "TriangleCount")
        GL.Uniform1(uTriangleCount, Me.triangleCount)

        'Dim CameraProjection As Integer = GL.GetUniformLocation(program, "CameraProjection")
        'GL.UniformMatrix4(CameraProjection, False, mtxProjection)
        'Dim CameraView As Integer = GL.GetUniformLocation(program, "CameraView")
        'GL.UniformMatrix4(CameraView, False, mtxView)

        Dim CameraPosition As Integer = GL.GetUniformLocation(program, "CameraPosition")
        GL.Uniform4(CameraPosition, cameraPos.X, cameraPos.Y, cameraPos.Z, cameraPos.W)

        Dim uSampler0 As Integer = GL.GetUniformLocation(program, "Sampler0")
        GL.Uniform1(uSampler0, 0) REM 第0个采样器
    End Sub

    Private bNormal As Boolean
    ''' <summary>
    ''' 递归遍历场景中的节点，把每个三角形添加到triangles列表中
    ''' </summary>
    ''' <param name="node">遞歸的父節點</param>
    ''' <param name="parentMtx">父變換矩陣</param>
    Public Sub RecursiveNode(node As Assimp.Node, parentMtx As Assimp.Matrix4x4)
        Dim program As Integer = ShaderManager.Active("RayTracing")
        Dim finalMtx As Assimp.Matrix4x4 = node.Transform * parentMtx

        Dim lightnames As List(Of String) = (From l In scene.Lights Select l.Name).ToList()
        If lightnames.Contains(node.Name) Then REM 节点是灯光
            Dim lightindex As Integer = lightnames.IndexOf(node.Name)
            Dim assimplight As Assimp.Light = scene.Lights(lightindex)
            Dim light As New Light(program, assimplight.Name)
            light.InitUniform()
            light.LightAttenuation = New Vector3(assimplight.AttenuationConstant,
                                                 assimplight.AttenuationLinear,
                                                 assimplight.AttenuationQuadratic)
            light.LightColor = ToVector3(assimplight.ColorDiffuse)
            light.LightDirection = ToVector3(New Matrix3x3(finalMtx) * assimplight.Direction)
            light.LightPosition = New Vector4(ToVector3(finalMtx * assimplight.Position), 1.0)
            light.SetupUniform()
        End If

        Dim cameranames As List(Of String) = (From c In scene.Cameras Select c.Name).ToList
        If cameranames.Contains(node.Name) Then REM 节点是摄像机
            Dim cameraindex As Integer = cameranames.IndexOf(node.Name)
            Dim cam As Assimp.Camera = scene.Cameras(cameraindex)
            Dim pos As Assimp.Vector3D = finalMtx * cam.Position
            Dim dir As Assimp.Vector3D = New Matrix3x3(finalMtx) * cam.Direction
            Dim up As Assimp.Vector3D = New Matrix3x3(finalMtx) * cam.Up
            cameraPos = New Vector4(ToVector3(pos), 1.0)
            mtxView = Matrix4.LookAt(cameraPos.Xyz, ToVector3(dir), ToVector3(up))
            mtxProjection = Matrix4.CreatePerspectiveFieldOfView(45 * Math.PI / 180, cam.AspectRatio, cam.ClipPlaneNear, cam.ClipPlaneFar)
            mtxModel = Matrix4.Identity
            HUDSurface.Singleton.MakeHUD(45, cam.AspectRatio, mtxView)
            HUDSurface.Singleton.SetupHUD()
        End If

        If node.HasMeshes Then
            For Each mi In node.MeshIndices
                Dim mesh As Assimp.Mesh = scene.Meshes(mi)
                REM 设置三角形跟网格映射，对应网格的变换 材质
                meshesInfo.Add(New MeshInfo(ToMatrix4(finalMtx), mesh.MaterialIndex))
                For fi = 0 To mesh.FaceCount - 1

                    Dim face As Assimp.Face = mesh.Faces(fi)
                    If face.IndexCount = 3 Then REM 仅限三角形
                        Dim A As Integer = face.Indices(0)
                        Dim B As Integer = face.Indices(1)
                        Dim C As Integer = face.Indices(2)

                        Dim tmp As Matrix4x4 = finalMtx
                        tmp.Inverse()
                        tmp.Transpose()
                        Dim finalNormalMtx As Matrix3x3 = New Matrix3x3(tmp) REM 法线矩阵 去掉模型矩阵的平移 缩放

                        Dim p0 As Assimp.Vector3D = finalMtx * mesh.Vertices(A)
                        Dim p1 As Assimp.Vector3D = finalMtx * mesh.Vertices(B)
                        Dim p2 As Assimp.Vector3D = finalMtx * mesh.Vertices(C)
                        REM 4float 32byte对齐

#Region "三角形信息"
                        REM p0,點的w分量非0
                        sceneData.WritePositionA(triangleCount, p0.X, p0.Y, p0.Z, 1)
                        sceneData.WritePositionB(triangleCount, p1.X, p1.Y, p1.Z, 1)
                        sceneData.WritePositionC(triangleCount, p2.X, p2.Y, p2.Z, 1)

                        Dim n0 As Assimp.Vector3D = finalNormalMtx * mesh.Normals(A)
                        Dim n1 As Assimp.Vector3D = finalNormalMtx * mesh.Normals(B)
                        Dim n2 As Assimp.Vector3D = finalNormalMtx * mesh.Normals(C)
                        n0.Normalize() : n1.Normalize() : n2.Normalize()

                        If bNormal Then
                            REM 法線 方向的w分量為0
                            sceneData.WriteNormalA(triangleCount, n0.X, n0.Y, n0.Z, 0)
                            sceneData.WriteNormalB(triangleCount, n1.X, n1.Y, n1.Z, 0)
                            sceneData.WriteNormalC(triangleCount, n2.X, n2.Y, n2.Z, 0)

                        Else
                            sceneData.WriteNormalA(triangleCount, 0, 0, 0, 0)
                            sceneData.WriteNormalB(triangleCount, 0, 0, 0, 0)
                            sceneData.WriteNormalC(triangleCount, 0, 0, 0, 0)
                        End If

                        If mesh.HasTextureCoords(0) Then
                            Dim t0 As Assimp.Vector3D = mesh.TextureCoordinateChannels(0)(A)
                            Dim t1 As Assimp.Vector3D = mesh.TextureCoordinateChannels(0)(B)
                            Dim t2 As Assimp.Vector3D = mesh.TextureCoordinateChannels(0)(C)
                            sceneData.WriteTexcoordA(triangleCount, t0.X, t0.Y)
                            sceneData.WriteTexcoordB(triangleCount, t1.X, t1.Y)
                            sceneData.WriteTexcoordC(triangleCount, t2.X, t2.Y)

                        Else
                            sceneData.WriteTexcoordA(triangleCount, 0, 0)
                            sceneData.WriteTexcoordB(triangleCount, 0, 0)
                            sceneData.WriteTexcoordC(triangleCount, 0, 0)
                        End If

                        REM 通过网格ID可以找到网格的变换矩阵，及其材质
                        'tData.Add(meshesInfo.Count - 1) REM 确定该三角形是属于哪一个mesh的 meshIndex
                        REM 虽然通过网格可以找到材质，但还是可以使用与材质不同的纹理
                        'Dim glTexID As Integer = textures.IndexOf(mesh.MaterialIndex)
                        'tData.Add(glTexID) REM 确定该三角形使用那个材质的贴图textureIndex

                        Dim mtl = scene.Materials(mesh.MaterialIndex)
                        sceneData.WriteInfo(triangleCount, meshesInfo.Count - 1,
                        mesh.MaterialIndex, 0)
#End Region
                        triangleCount += 1
                    End If
                Next
            Next
        End If

        If node.HasChildren Then
            For Each child As Assimp.Node In node.Children
                RecursiveNode(child, finalMtx)
            Next
        End If
    End Sub

    Public Function GetTriangleCount() As Integer
        Dim r As Integer = 0
        For Each m In scene.Meshes
            r += m.FaceCount
        Next
        Return r
    End Function

    Public vbo As Integer
    Public bufferOutput As Integer
    Public Sub CreateVBO(data() As Single)
        vbo = GL.GenBuffer()
        GL.BindBuffer(BufferTarget.ArrayBuffer, vbo)
        GL.BufferData(BufferTarget.ArrayBuffer, data.Length * 4, data, BufferUsageHint.StaticDraw)
        GL.BindBuffer(BufferTarget.ArrayBuffer, 0)

        bufferOutput = GL.GenBuffer()
        GL.BindBuffer(BufferTarget.ArrayBuffer, bufferOutput)
        GL.BufferData(BufferTarget.ArrayBuffer, data.Length * 4, data, BufferUsageHint.StaticDraw)
        GL.BindBuffer(BufferTarget.ArrayBuffer, 0)

    End Sub

    Private vao As Integer
    Public Sub CreateVAO()
        Dim id As Integer = ShaderManager.Active("GBuffer")
        vao = GL.GenVertexArray()
        GL.BindVertexArray(vao)
        GL.BindBuffer(BufferTarget.ArrayBuffer, vbo)
        GL.VertexAttribPointer(0, 3, VertexAttribPointerType.Float, False, triangleCount, 0)
        GL.VertexAttribPointer(1, 3, VertexAttribPointerType.Float, False, triangleCount, triangleCount * 3)
        GL.VertexAttribPointer(2, 2, VertexAttribPointerType.Float, False, triangleCount, triangleCount * 6)
        GL.BindVertexArray(0)
    End Sub

    Public meshesInfo As New List(Of MeshInfo)
    ''' <summary>
    ''' 把模型中的所有三角形存到一張圖片轉成一張紋理
    ''' </summary>
    Public Function BuildTrianglesDatas() As Single()
        Dim finalMtx As Assimp.Matrix4x4 = Matrix4x4.Identity
        REM 重置数据
        triangleCount = 0
        meshesInfo.Clear()
        bNormal = True
        RecursiveNode(scene.RootNode, finalMtx)

        'REM 一行1024个纹素 32 个三角形 4096个float
        REM 一行8个纹素 1个三角形 32个float
        Dim floatbuf As Single() = sceneData.GetData()
        CreateVBO(floatbuf)
        'SetupVerticesTexture(floatbuf)
        UpdateVerticeTextureBuffer()
        SetupMeshesInfo()
        Return floatbuf
    End Function

    Public Sub SetupMeshesInfo()
        Dim program As Integer = ShaderManager.Active("RayTracing")
        For i = 0 To meshesInfo.Count - 1
            Dim infoTransform As Integer = GL.GetUniformLocation(program, "meshesInfo[" + i.ToString() + "].LocalTransform")
            GL.UniformMatrix4(infoTransform, False, meshesInfo(i).localTransform) REM 三角形对应的网格要进行的变换
            Dim infoMaterialID As Integer = GL.GetUniformLocation(program, "meshesInfo[" + i.ToString() + "].MaterialID")
            GL.Uniform1(infoMaterialID, meshesInfo(i).materialID) REM 三角形对应的网格的纹理材质
        Next
    End Sub

    Public Sub SetupVerticesTexture(varr() As Single)
        GL.ActiveTexture(TextureUnit.Texture0)
        GL.BindTexture(TextureTarget.Texture2D, dataTexture)
        GL.TexParameter(TextureTarget.Texture2D, TextureParameterName.TextureMagFilter, TextureMagFilter.Nearest)
        GL.TexParameter(TextureTarget.Texture2D, TextureParameterName.TextureMinFilter, TextureMinFilter.Nearest)

        ' 一行一个三角形的局部信息，共10行，这样的纹理最大存16384个
        GL.TexImage2D(TextureTarget.Texture2D, 0, PixelInternalFormat.Rgba32f, triangleCount, 10, 0, PixelFormat.Rgba, PixelType.Float, varr)
        GL.TexBuffer(TextureBufferTarget.TextureBuffer, SizedInternalFormat.Rgba32f, bufferOutput)
    End Sub

    Public Sub UpdateVerticeTextureBuffer()
        GL.ActiveTexture(TextureUnit.Texture0)
        GL.BindTexture(TextureTarget.TextureBuffer, dataTexture)
        GL.TexBuffer(TextureBufferTarget.TextureBuffer, SizedInternalFormat.Rgba32f, bufferOutput)
    End Sub

    Public Sub SetupTextureHandle(program As Integer)
        If BindlessTextureSupported() Then
            Dim ubindless = GL.GetUniformLocation(program, "bindless")
            Dim hsize = RayTracing.sceneManager.texturehandles.Length
            GL.Arb.UniformHandle(ubindless, hsize, RayTracing.sceneManager.texturehandles)
        End If
    End Sub

    Public Sub SetupMaterial()
        Dim program As Integer = ShaderManager.Active("RayTracing")
        For i = 0 To scene.Materials.Count - 1
            Dim mtl As New Material(program, scene.Materials(i).Name)
            Dim v = scene.Materials(i).GetAllProperties()

            mtl.InitUniform()
            mtl.MaterialDiffuse = ToVector4(scene.Materials(i).ColorDiffuse)
            mtl.MaterialEmissive = ToVector4(scene.Materials(i).ColorEmissive)
            mtl.MaterialSpecular = ToVector4(scene.Materials(i).ColorSpecular)
            mtl.MaterialOpacity = scene.Materials(i).Opacity
            mtl.MaterialSpecularShiness = scene.Materials(i).Shininess
            mtl.MaterialReflectivity = scene.Materials(i).Reflectivity
            mtl.SetupUniform()
        Next
    End Sub

    Public Function ReadBuffer(buf As Integer) As Single()
        GL.BindBuffer(BufferTarget.CopyReadBuffer, buf)
        Dim bufsize As Integer
        GL.GetBufferParameter(BufferTarget.CopyReadBuffer, BufferParameterName.BufferSize, bufsize)
        Dim data() = New Single(bufsize / 4 - 1) {}
        GL.GetBufferSubData(Of Single)(BufferTarget.CopyReadBuffer, IntPtr.Zero, bufsize, data)
        Return data
    End Function

    Public Shared Function GetBoxAABB(data() As Single) As Boundingbox
        Dim result As Boundingbox
        result.Type = BoxType.Aabb
        Dim minX, minY, minZ, maxX, maxY, maxZ As Single
        For i = 0 To data.Length - 1 Step 4
            Dim x = data(i + 0)
            Dim y = data(i + 1)
            Dim z = data(i + 2)
            If x > maxX Then maxX = x
            If y > maxY Then maxY = y
            If z > maxZ Then maxZ = z
            If x < minX Then minX = x
            If y < minY Then minY = y
            If z < minZ Then minZ = z
        Next
        result.Max = New Vector3(maxX, maxY, maxZ)
        result.Min = New Vector3(minX, minY, minZ)
        Return result
    End Function

    Public Shared Function GetBoxSphere(data() As Single) As Boundingbox
        Dim result As Boundingbox
        result.Type = BoxType.Sphere
        Dim r As Single = 0
        Dim c As Vector3 = New Vector3(data(0), data(1), data(2))
        Dim num = data.Length \ 4
        For i = 4 To data.Length - 1 Step 4
            Dim x = data(i + 0)
            Dim y = data(i + 1)
            Dim z = data(i + 2)
            Dim index = i \ 4
            Dim p As New Vector3(x, y, z)
            Dim dist = (c - p).LengthFast
            If dist < r Then Continue For

            r = (r + dist) / 2
            c = p + Vector3.Normalize(c - p) * r
        Next
        result.Center = c
        result.Radius = r
        Return result
    End Function

    '方差 = (X[0] - X[avg])^2 + ... + (X[n] - X[avg])^2) / (n - 1)
    '标准差 = sqrt(方差)
    '协方差 Cov(x, y) = Cov(y, x) = sum((X[n] - X[avg]) * (Y[n] - Y[avg])) / (n - 1)

    'det
    '[a, b] = ad - bc
    '[c, d]

    '<=>
    '[a, b, 0]
    '[c, d, 0] = a(d*1 - 0*0) - b(c*1- 0*0) + 0(c*0- d*0)
    '[0, 0, 1] = ad - bc

    '[a, b, c] = a(ei-hf) - b(di-gf) +c(dh-eg)
    '[d, e, f]
    '[g, h, i]

    '[a,b,c,d] = a*[f(kp-lo)-g(jp-ln)+h(jo-kn)] - b[e(kp-lo)-g(ip-lm)+h(io-km)] +
    '[e,f,g,h]   c*[e(jp-ln)-f(ip-lm)+h(in-jm)] - d[e(jo-kn)-f(io-km)+g(in-jm)]
    '[i,j,k,l]
    '[m,n,o,p]

    Public Shared Function GetBoxObb(data() As Single) As Boundingbox
        Dim result As Boundingbox
        result.Type = BoxType.Obb
        REM  计算平均值
        Dim avg = Vector3.Zero
        For i = 0 To data.Length - 1 Step 4
            Dim x = data(i + 0)
            Dim y = data(i + 1)
            Dim z = data(i + 2)
            Dim p As New Vector3(x, y, z)
            avg += p
        Next
        Dim n As Integer = data.Length / 4
        avg = avg / n

        REM PCA 主成分分析 https://blog.csdn.net/program_developer/article/details/80632779
        REM obb 生成 https://blog.csdn.net/qing101hua/article/details/53100112
        REM 计算协方差，协方差Cov(x,y)可以描述一个n维数据集在其中某两个维度的相关性
        REM 协方差矩阵 mtxCov = 
        REM Cov(x,x) Cov(x,y) Cov(x,z)
        REM Cov(y,x) Cov(y,y) Cov(y,z)
        REM Cov(z,x) Cov(z,y) Cov(z,z)
        REM 其中 Cov(x,y) = Cov(y,x) Cov(x,z) = Cov(z,x) Cov(y,z) = Cov(z,y)
        REM 有偏估计 1 / n  sum(Xi - avg)
        REM 无偏估计 n / (n-1) * 有偏估计 =>  1 / (n-1) SUM(Xi - avg)
        REM 协方差  1 / (n-1) sum((Xi - avg.x)(Yi - avg.y))
        Dim covXX, covXY, covXZ, covYY, covYZ, covZZ As Single
        For i = 0 To data.Length - 1 Step 4
            Dim x = data(i + 0)
            Dim y = data(i + 1)
            Dim z = data(i + 2)
            covXX += (x - avg.X) * (x - avg.X)
            covXY += (x - avg.X) * (y - avg.Y)
            covXZ += (x - avg.X) * (z - avg.Z)
            covYY += (y - avg.Y) * (y - avg.Y)
            covYZ += (y - avg.Y) * (z - avg.Z)
            covZZ += (z - avg.Z) * (z - avg.Z)
        Next
        covXX /= (n - 1)
        covXY /= (n - 1)
        covXZ /= (n - 1)
        covYY /= (n - 1)
        covYZ /= (n - 1)
        covZZ /= (n - 1)
        Dim matCov As New Matrix3(covXX, covXY, covXZ, covXY, covYY, covYZ, covXZ, covYZ, covZZ)
        'http://www.360doc.com/content/17/0727/13/28378250_674516083.shtml jocobianIteration to solve eigen vector
        Dim evalue = EigenValue(matCov)
        Dim evector = EigenVector(matCov, evalue)

        Array.Sort(evalue, evector, New SortAbsDesc())
        Dim primeAxis As Vector3 = evector(0)

        Dim secondAxis As Vector3
        Dim vAvg = New Vector3(avg.X, avg.Y, avg.Z)
        If evalue.Length > 1 Then
            secondAxis = evector(1)
        Else
            'test secondAxis
            Dim cov1XX, cov1XY, cov1XZ, cov1YY, cov1YZ, cov1ZZ As Single
            For i = 0 To data.Length - 1 Step 4
                Dim x = data(i + 0)
                Dim y = data(i + 1)
                Dim z = data(i + 2)
                Dim tmp = New Vector3(x, y, z) - vAvg
                Dim v = tmp - primeAxis * Vector3.Dot(primeAxis, tmp)
                cov1XX += v.X * v.X
                cov1XY += v.X * v.Y
                cov1XZ += v.X * v.Z
                cov1YY += v.Y * v.Y
                cov1YZ += v.Y * v.Z
                cov1ZZ += v.Z * v.Z
            Next
            cov1XX /= (n - 1)
            cov1XY /= (n - 1)
            cov1XZ /= (n - 1)
            cov1YY /= (n - 1)
            cov1YZ /= (n - 1)
            cov1ZZ /= (n - 1)

            Dim matCov1 As New Matrix3(cov1XX, cov1XY, cov1XZ, cov1XY, cov1YY, cov1YZ, cov1XZ, cov1YZ, cov1ZZ)
            Dim evalue1 = EigenValue(matCov1)
            Dim evector1 = EigenVector(matCov1, evalue1)
            Array.Sort(evalue1, evector1, New SortAbsDesc())
            secondAxis = evector1(0)
        End If

        Dim thirdAxis = Vector3.Cross(primeAxis, secondAxis).Normalized()
        Dim minX, minY, minZ, maxX, maxY, maxZ As Single
        For i = 0 To data.Length - 1 Step 4
            Dim x = data(i + 0)
            Dim y = data(i + 1)
            Dim z = data(i + 2)
            Dim delta = New Vector3(x, y, z) - vAvg
            Dim dx = Vector3.Dot(delta, primeAxis)
            Dim dy = Vector3.Dot(delta, secondAxis)
            Dim dz = Vector3.Dot(delta, thirdAxis)
            If dx > maxX Then maxX = dx
            If dy > maxY Then maxY = dy
            If dz > maxZ Then maxZ = dz
            If dx < minX Then minX = dx
            If dy < minY Then minY = dy
            If dz < minZ Then minZ = dz
        Next
        result.Center = vAvg
        result.Min = New Vector3(minX, minY, minZ)
        result.Max = New Vector3(maxX, maxY, maxZ)

        Dim aa = Helper.GetRotationOfCoordinateSystem(primeAxis, secondAxis, thirdAxis)
        result.Rotation = OpenTK.Quaternion.FromAxisAngle(aa.Xyz, aa.W)

        result.vPrime = result.Center + Vector3.UnitX * 8
        result.vSecond = result.Center + Vector3.UnitY * 4
        result.vThird = result.Center + Vector3.UnitZ * 2
        Return result
    End Function

End Class

Public Enum BoxType
    Sphere
    Aabb
    Obb
End Enum

Public Structure Boundingbox
    Public Type As BoxType
    Public Center As Vector3
    Public Radius As Single
    Public Min As Vector3
    Public Max As Vector3
    Public Rotation As OpenTK.Quaternion
    Public vPrime As Vector3
    Public vSecond As Vector3
    Public vThird As Vector3
End Structure

Public Class ConvexHull
    Public edgeVertices As List(Of Vector3)

    Public Sub New(vts As Vector3())
        Dim maxIX As Integer = 0
        Dim maxIY As Integer = 0
        Dim maxIZ As Integer = 0
        Dim minIX As Integer = 0
        Dim minIY As Integer = 0
        Dim minIZ As Integer = 0
        For i = 0 To vts.Length - 1
            If vts(i).X > vts(maxIX).X Then
                maxIX = i
            End If
            If vts(i).Y > vts(maxIY).Y Then
                maxIY = i
            End If
            If vts(i).Z > vts(maxIZ).Z Then
                maxIZ = i
            End If
            If vts(i).X < vts(minIX).X Then
                minIX = i
            End If
            If vts(i).Y < vts(minIY).Y Then
                minIY = i
            End If
            If vts(i).Z < vts(minIZ).Z Then
                minIZ = i
            End If
        Next

    End Sub
End Class

Public Interface IWriteSceneData
    Sub WritePositionA(index As Integer, x As Single, y As Single, z As Single, w As Single)
    Sub WritePositionB(index As Integer, x As Single, y As Single, z As Single, w As Single)
    Sub WritePositionC(index As Integer, x As Single, y As Single, z As Single, w As Single)
    Sub WriteNormalA(index As Integer, x As Single, y As Single, z As Single, w As Single)
    Sub WriteNormalB(index As Integer, x As Single, y As Single, z As Single, w As Single)
    Sub WriteNormalC(index As Integer, x As Single, y As Single, z As Single, w As Single)
    Sub WriteTexcoordA(index As Integer, x As Single, y As Single, Optional z As Single = 0, Optional w As Single = 0)
    Sub WriteTexcoordB(index As Integer, x As Single, y As Single, Optional z As Single = 0, Optional w As Single = 0)
    Sub WriteTexcoordC(index As Integer, x As Single, y As Single, Optional z As Single = 0, Optional w As Single = 0)
    Sub WriteInfo(index As Integer, meshindex As Integer, mtlindex As Integer, tagindex As Integer)
    Function GetData() As Array
End Interface

Public Class InterleaveSceneData
    Implements IWriteSceneData
    Public data() As Single

    Public Sub New()
        data = New Single(SceneManager.MaxTriangles * 40 - 1) {}
    End Sub

    Public Sub New(size As Integer)
        data = New Single(size * 40 - 1) {}
    End Sub

    Public Sub WritePositionA(index As Integer, x As Single, y As Single, z As Single, w As Single) Implements IWriteSceneData.WritePositionA
        data(index * 40 + 0) = x
        data(index * 40 + 1) = y
        data(index * 40 + 2) = z
        data(index * 40 + 3) = w
    End Sub

    Public Sub WritePositionB(index As Integer, x As Single, y As Single, z As Single, w As Single) Implements IWriteSceneData.WritePositionB
        data(index * 40 + 4 + 0) = x
        data(index * 40 + 4 + 1) = y
        data(index * 40 + 4 + 2) = z
        data(index * 40 + 4 + 3) = w
    End Sub

    Public Sub WritePositionC(index As Integer, x As Single, y As Single, z As Single, w As Single) Implements IWriteSceneData.WritePositionC
        data(index * 40 + 8 + 0) = x
        data(index * 40 + 8 + 1) = y
        data(index * 40 + 8 + 2) = z
        data(index * 40 + 8 + 3) = w
    End Sub

    Public Sub WriteNormalA(index As Integer, x As Single, y As Single, z As Single, w As Single) Implements IWriteSceneData.WriteNormalA
        data(index * 40 + 12 + 0) = x
        data(index * 40 + 12 + 1) = y
        data(index * 40 + 12 + 2) = z
        data(index * 40 + 12 + 3) = w
    End Sub

    Public Sub WriteNormalB(index As Integer, x As Single, y As Single, z As Single, w As Single) Implements IWriteSceneData.WriteNormalB
        data(index * 40 + 16 + 0) = x
        data(index * 40 + 16 + 1) = y
        data(index * 40 + 16 + 2) = z
        data(index * 40 + 16 + 3) = w
    End Sub

    Public Sub WriteNormalC(index As Integer, x As Single, y As Single, z As Single, w As Single) Implements IWriteSceneData.WriteNormalC
        data(index * 40 + 20 + 0) = x
        data(index * 40 + 20 + 1) = y
        data(index * 40 + 20 + 2) = z
        data(index * 40 + 20 + 3) = w
    End Sub

    Public Sub WriteTexcoordA(index As Integer, x As Single, y As Single, Optional z As Single = 0, Optional w As Single = 0) Implements IWriteSceneData.WriteTexcoordA
        data(index * 40 + 24 + 0) = x
        data(index * 40 + 24 + 1) = y
        data(index * 40 + 24 + 2) = z
        data(index * 40 + 24 + 3) = w
    End Sub

    Public Sub WriteTexcoordB(index As Integer, x As Single, y As Single, Optional z As Single = 0, Optional w As Single = 0) Implements IWriteSceneData.WriteTexcoordB
        data(index * 40 + 28 + 0) = x
        data(index * 40 + 28 + 1) = y
        data(index * 40 + 28 + 2) = z
        data(index * 40 + 28 + 3) = w
    End Sub

    Public Sub WriteTexcoordC(index As Integer, x As Single, y As Single, Optional z As Single = 0, Optional w As Single = 0) Implements IWriteSceneData.WriteTexcoordC
        data(index * 40 + 32 + 0) = x
        data(index * 40 + 32 + 1) = y
        data(index * 40 + 32 + 2) = z
        data(index * 40 + 32 + 3) = w
    End Sub

    Public Sub WriteMaterailIndex(index As Integer, meshindex As Integer, mtlindex As Integer, tagindex As Integer) Implements IWriteSceneData.WriteInfo
        data(index * 40 + 36) = index
        data(index * 40 + 37) = meshindex
        data(index * 40 + 38) = mtlindex
        data(index * 40 + 39) = tagindex
    End Sub

    Public Function GetData() As Array Implements IWriteSceneData.GetData
        Return data
    End Function
End Class

Public Class SeparateSceneData
    Implements IWriteSceneData

    Public positionsA() As Single
    Public positionsB() As Single
    Public positionsC() As Single

    Public normalsA() As Single
    Public normalsB() As Single
    Public normalsC() As Single

    Public texcoordsA() As Single
    Public texcoordsB() As Single
    Public texcoordsC() As Single

    Public Info() As Single

    Private triangleCount As Integer
    Public Sub New()
        triangleCount = SceneManager.MaxTriangles
        ReDim positionsA(SceneManager.MaxTriangles * 4 - 1)
        ReDim positionsB(SceneManager.MaxTriangles * 4 - 1)
        ReDim positionsC(SceneManager.MaxTriangles * 4 - 1)

        ReDim normalsA(SceneManager.MaxTriangles * 4 - 1)
        ReDim normalsB(SceneManager.MaxTriangles * 4 - 1)
        ReDim normalsC(SceneManager.MaxTriangles * 4 - 1)

        ReDim texcoordsA(SceneManager.MaxTriangles * 4 - 1)
        ReDim texcoordsB(SceneManager.MaxTriangles * 4 - 1)
        ReDim texcoordsC(SceneManager.MaxTriangles * 4 - 1)
        ReDim Info(SceneManager.MaxTriangles * 4 - 1)
    End Sub

    Public Sub New(size As Integer)
        triangleCount = size
        ReDim positionsA(size * 4 - 1)
        ReDim positionsB(size * 4 - 1)
        ReDim positionsC(size * 4 - 1)

        ReDim normalsA(size * 4 - 1)
        ReDim normalsB(size * 4 - 1)
        ReDim normalsC(size * 4 - 1)

        ReDim texcoordsA(size * 4 - 1)
        ReDim texcoordsB(size * 4 - 1)
        ReDim texcoordsC(size * 4 - 1)
        ReDim Info(size * 4 - 1)
    End Sub

    Public Sub WritePositionA(index As Integer, x As Single, y As Single, z As Single, w As Single) Implements IWriteSceneData.WritePositionA
        positionsA(index * 4 + 0) = x
        positionsA(index * 4 + 1) = y
        positionsA(index * 4 + 2) = z
        positionsA(index * 4 + 3) = w
    End Sub

    Public Sub WritePositionB(index As Integer, x As Single, y As Single, z As Single, w As Single) Implements IWriteSceneData.WritePositionB
        positionsB(index * 4 + 0) = x
        positionsB(index * 4 + 1) = y
        positionsB(index * 4 + 2) = z
        positionsB(index * 4 + 3) = w
    End Sub

    Public Sub WritePositionC(index As Integer, x As Single, y As Single, z As Single, w As Single) Implements IWriteSceneData.WritePositionC
        positionsC(index * 4 + 0) = x
        positionsC(index * 4 + 1) = y
        positionsC(index * 4 + 2) = z
        positionsC(index * 4 + 3) = w
    End Sub

    Public Sub WriteNormalA(index As Integer, x As Single, y As Single, z As Single, w As Single) Implements IWriteSceneData.WriteNormalA
        normalsA(index * 4 + 0) = x
        normalsA(index * 4 + 1) = y
        normalsA(index * 4 + 2) = z
        normalsA(index * 4 + 3) = w
    End Sub

    Public Sub WriteNormalB(index As Integer, x As Single, y As Single, z As Single, w As Single) Implements IWriteSceneData.WriteNormalB
        normalsB(index * 4 + 0) = x
        normalsB(index * 4 + 1) = y
        normalsB(index * 4 + 2) = z
        normalsB(index * 4 + 3) = w
    End Sub

    Public Sub WriteNormalC(index As Integer, x As Single, y As Single, z As Single, w As Single) Implements IWriteSceneData.WriteNormalC
        normalsC(index * 4 + 0) = x
        normalsC(index * 4 + 1) = y
        normalsC(index * 4 + 2) = z
        normalsC(index * 4 + 3) = w
    End Sub

    Public Sub WriteTexcoordA(index As Integer, x As Single, y As Single, Optional z As Single = 0, Optional w As Single = 0) Implements IWriteSceneData.WriteTexcoordA
        texcoordsA(index * 4 + 0) = x
        texcoordsA(index * 4 + 1) = y
        texcoordsA(index * 4 + 2) = z
        texcoordsA(index * 4 + 3) = w
    End Sub

    Public Sub WriteTexcoordB(index As Integer, x As Single, y As Single, Optional z As Single = 0, Optional w As Single = 0) Implements IWriteSceneData.WriteTexcoordB
        texcoordsB(index * 4 + 0) = x
        texcoordsB(index * 4 + 1) = y
        texcoordsB(index * 4 + 2) = z
        texcoordsB(index * 4 + 3) = w
    End Sub

    Public Sub WriteTexcoordC(index As Integer, x As Single, y As Single, Optional z As Single = 0, Optional w As Single = 0) Implements IWriteSceneData.WriteTexcoordC
        texcoordsC(index * 4 + 0) = x
        texcoordsC(index * 4 + 1) = y
        texcoordsC(index * 4 + 2) = z
        texcoordsC(index * 4 + 3) = w
    End Sub

    Public Sub WriteInfo(index As Integer, meshindex As Integer, mtlindex As Integer, tagIndex As Integer) Implements IWriteSceneData.WriteInfo
        Info(index * 4 + 0) = index
        Info(index * 4 + 1) = meshindex
        Info(index * 4 + 2) = mtlindex
        Info(index * 4 + 3) = tagIndex
    End Sub

    Private tmpData As Single()
    Public Function GetData() As Array Implements IWriteSceneData.GetData
        If tmpData IsNot Nothing Then Return tmpData
        Dim data = New Single(triangleCount * 40 - 1) {}
        Array.Copy(positionsA, 0, data, 0 * triangleCount * 4, triangleCount * 4)
        Array.Copy(positionsB, 0, data, 1 * triangleCount * 4, triangleCount * 4)
        Array.Copy(positionsC, 0, data, 2 * triangleCount * 4, triangleCount * 4)
        Array.Copy(normalsA, 0, data, 3 * triangleCount * 4, triangleCount * 4)
        Array.Copy(normalsB, 0, data, 4 * triangleCount * 4, triangleCount * 4)
        Array.Copy(normalsC, 0, data, 5 * triangleCount * 4, triangleCount * 4)
        Array.Copy(texcoordsA, 0, data, 6 * triangleCount * 4, triangleCount * 4)
        Array.Copy(texcoordsB, 0, data, 7 * triangleCount * 4, triangleCount * 4)
        Array.Copy(texcoordsC, 0, data, 8 * triangleCount * 4, triangleCount * 4)
        Array.Copy(Info, 0, data, 9 * triangleCount * 4, triangleCount * 4)
        tmpData = data
        Return data
    End Function

    Private tmpVertices As Single()
    Public Function GetVertices() As Array
        If tmpVertices IsNot Nothing Then Return tmpVertices
        Dim data = New Single(triangleCount * 36 - 1) {}
        Array.Copy(positionsA, 0, data, 0 * triangleCount * 4, triangleCount * 4)
        Array.Copy(positionsB, 0, data, 1 * triangleCount * 4, triangleCount * 4)
        Array.Copy(positionsC, 0, data, 2 * triangleCount * 4, triangleCount * 4)
        Array.Copy(normalsA, 0, data, 3 * triangleCount * 4, triangleCount * 4)
        Array.Copy(normalsB, 0, data, 4 * triangleCount * 4, triangleCount * 4)
        Array.Copy(normalsC, 0, data, 5 * triangleCount * 4, triangleCount * 4)
        Array.Copy(texcoordsA, 0, data, 6 * triangleCount * 4, triangleCount * 4)
        Array.Copy(texcoordsB, 0, data, 7 * triangleCount * 4, triangleCount * 4)
        Array.Copy(texcoordsC, 0, data, 8 * triangleCount * 4, triangleCount * 4)
        tmpVertices = data
        Return data
    End Function

    Public Function IndexOfTriangle(index As Integer) As TriangleData
        Dim r As TriangleData
        r.p1 = New Vector4(positionsA(index * 4 + 0), positionsA(index * 4 + 1), positionsA(index * 4 + 2), positionsA(index * 4 + 3))
        r.p2 = New Vector4(positionsB(index * 4 + 0), positionsB(index * 4 + 1), positionsB(index * 4 + 2), positionsB(index * 4 + 3))
        r.p3 = New Vector4(positionsC(index * 4 + 0), positionsC(index * 4 + 1), positionsC(index * 4 + 2), positionsC(index * 4 + 3))
        r.n1 = New Vector4(normalsA(index * 4 + 0), normalsA(index * 4 + 1), normalsA(index * 4 + 2), normalsA(index * 4 + 3))
        r.n2 = New Vector4(normalsB(index * 4 + 0), normalsB(index * 4 + 1), normalsB(index * 4 + 2), normalsB(index * 4 + 3))
        r.n3 = New Vector4(normalsC(index * 4 + 0), normalsC(index * 4 + 1), normalsC(index * 4 + 2), normalsC(index * 4 + 3))
        r.t1 = New Vector4(texcoordsA(index * 4 + 0), texcoordsA(index * 4 + 1), texcoordsA(index * 4 + 2), texcoordsA(index * 4 + 3))
        r.t2 = New Vector4(texcoordsB(index * 4 + 0), texcoordsB(index * 4 + 1), texcoordsB(index * 4 + 2), texcoordsB(index * 4 + 3))
        r.t3 = New Vector4(texcoordsC(index * 4 + 0), texcoordsC(index * 4 + 1), texcoordsC(index * 4 + 2), texcoordsC(index * 4 + 3))
        r.Info = New Vector4(Info(index * 4 + 0), Info(index * 4 + 1), Info(index * 4 + 2), Info(index * 4 + 3))
        Return r
    End Function

End Class

Public Structure TriangleData
    Public p1 As Vector4
    Public p2 As Vector4
    Public p3 As Vector4
    Public n1 As Vector4
    Public n2 As Vector4
    Public n3 As Vector4
    Public t1 As Vector4
    Public t2 As Vector4
    Public t3 As Vector4
    Public Info As Vector4
End Structure